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Traditionally, finite differences and finite element methods have been by many regarded as the 
basic tools for obtaining numerical solutions in a variety of quantum mechanical problems emerging 
in atomic, nuclear and particle physics, astrophysics, quantum chemistry, etc. In recent years, 
however, an alternative technique based on the semi-spectral methods has focused considerable 
attention. The purpose of this work is first to provide the necessary tools and subsequently examine 
the efficiency of this method in quantum mechanical applications. Restricting our interest to time 
independent two-body problems, we obtained the continuous and discrete spectrum solutions of the 
underlying Schrodinger or Lippmann-Schwinger equations in both, the coordinate and momentum 
, space. In all of the numerically studied examples we had no difficulty in achieving the machine 

■ accuracy and the semi-spectral method showed exponential convergence combined with excellent 

' numerical stability. 

(N 
(N 

I. INTRODUCTION 

^ ' The exact solutions of the Schrodinger equation are known to exist only in a few rather ideahzed examples. Since 
, in all remaining cases one has to be content with the numerically generated solutions, the importance of numerical 
' methods in quantum mechanics can hardly be overestimated. Broadly speaking, one may distinguish two schools of 
thought Jj,2]: (i) the local approach comprising the finite differences and finite element methods, and, (ii) the global 
approach encompassing a variety of spectral methods. The local methods usually require much larger grids than other 
methods as the convergence of these schemes is generally algebraic. By contrast, the global methods have been shown 
to provide exponential convergence delivering smooth solutions which in most cases are preferable. 
' The local methods approximate the unknown function by a sequence of overlapping low order polynomials which 
are used to interpolate the function at a small sub-set of the grid points. The derivative of the unknown function is 
then approximated by the derivative of the interpolating polynomial so that the former takes the form of a weighted 
, sum of the function values at the interpolation points. The method is designed to be exact for polynomials of low 
' order and the form of the difference quotient is chosen such that certain order of the truncation error is maintained 
^ I (in the lowest order du{x) /dx w [u{x + h) — u{x — h)]/ {2h) + 0{h^) where /i is a small grid spacing). In the simplest 
version this procedure converts a differential equation into a three-term recurrence, as may be exemplified by the 
^ popular Numerov algorithm easy to program and implement but relatively low in accuracy. 
• I— I , The strategy of the global or spectral methods is very different. Here, the unknown function is sought in the form 
of a generalized Fourier expansion, truncated for practical reasons, using global basis functions which usually are 
^ , polynomials or trigonometric polynomials of a high degree. This approach may be viewed as the limit of the finite 
- - ' difference method whose order of accuracy has been pushed to infinity. 

The semi-spectral Chebyshev method - a prominent representative of the global approach - provides effective means 
to solve a generic equation of the form L u{x) = f{x) where the independent variable x belongs to some interval [a, 5], 
L is a differential or integral operator, f{x) is a given function and u{x) is the unknown function satisfying some, as 
yet unspecified boundary conditions. For the sake of simplicity of the exposition we consider just a single independent 
variable x. The underlying idea is exceedingly simple: to obtain the solution, first the problem is discretized by 
introducing a grid consisting of N points ti £ [a, 6], z = 1, 2, . . . N , then the underlying equation is solved exactly at all 
of the grid points, and, finally with the solutions u{ti) in hand, an interpolating procedure is invoked to get u{x) at an 
arbitrarily chosen x value. Obviously, the key issue here would be the choice of the grid and the interpolation method 
but it turns out that the scheme based on Chebyshev polynomials provides the best offer on both accounts. It would 
be instructive to compare this approach with other methods and put it in a historical perspective but this is not so 
simple because the abundant mathematical literature is lacking a generally accepted standardization and quite 
often different methods bear the same tag, or else, the same thing appears under different names. Accordingly, the 
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semi-spectral method might be called orthogonal collocation, Petrov-Galerkin, or even Rayleigh-Ritz method. The 
common point of departure in all of the methods is the effort to minimize the error distribution given by the reminder 
or the residue function R{x) = Lu{x) — f{x) and there is a number of possible strategies to achieve this goal. It is 
convenient to seek the solution as a superposition of some trial functions (x) 



N 

u{x) =:^^Ci(f>,{x), (1) 



where Ci are hitherto unknown generalized Fourier expansion coefficients. The trial functions are expected to be easy 
to get by apart from the usual indispensable requirements of completeness and orthogonality. Local methods like 
finite difference or finite element methods use for this purpose low order partially overlapping polynomials defined 
on a grid whereas the global methods employ trigonometric functions (Fourier methods) or high order polynomials 
spanning over the whole domain. The global methods are also known as the spectral methods. In order to minimize 
the residue function another set of the so called test functions Xi{^) is needed together with a definition of a scalar 
product of two functions f{x) and g{x) which is taken to be given as 

(/,5)= f f{x)g{x)p{x)dx, (2) 



where p is a weighting function and the interval [a, b] specifies the domain. When u were the exact solution the residue 
function would identically vanish but with an approximate solution we can only try making the error distribution as 
small as possible by minimizing the residue function. This can be effected by stipulating the conditions (i?, Xi) — 

0, i ^ 1,2, . . . N . Indeed, since the test functions form a complete set, the vector R orthogonal to all members of such 
a set must be a null vector. The above orthogonality conditions lead to a linear set of algebraic equations 

N 

^(x,,L0.)c. = (x,,/), j = l,2,...iV (3) 

i=l 

for the expansion coefficients Ci which can be solved by standard methods. It will be convenient to distinguish three 
possible implementations of the above scheme. 

(i) Semi-spectral method. A grid ti is introduced and the test functions are adopted in the form Xi(.T) = S(x — ti). 

(ii) Galerkin, or spectral method. The test functions and the trial functions are identical Xi{^) ~ 4'i{x) fo^' 

(iii) Petrov-Galerkin method. The test functions and the trial fimctions are different. 

It has to be noted that the grid is indispensable only in the semi-spectral method whereas the remaining two methods 
can do without it. In the following we shall concentrate on the semi-spectral and the Galerkin method. 

As far as the basis functions are concerned, it is a common practice taking for 4>j{x) any of the sets of the orthogonal 
polynomials collected in j^l- Indeed, since all of them result as eigenf unctions of a certain Sturm-Liouville problem 
their completeness and orthogonality are guaranteed. There are also additional advantages following from taking 
polynomials as the expansion basis: firstly, this immediately sets an important correspondence linking the generalized 
Fourier expansion with an interpolating scheme, and, secondly, it simplifies all the necessary integrations by using the 
highly accurate Gauss integration scheme. We shall take on both of these issues, assuming that from now on (cc) 
stands for a polynomial of the order {j — 1). For an approximation of the order A'' and assigned polynomial type the 
grid points tj, j = 1,2, . . . N will be identical with the Gauss abscissas which are given as the zeros of the polynomial 
</>Af+i(a;). 

Integration. The Gauss integration formula, takes the form 

b N 
u{x) p{x) dx — Wj u(tj) (4) 

and is known to be exact if u{x) is a polynomial whose order is not greater than (2 A^ + 1). The weights Wj, j = 

1, 2, . . . A^ are obtained by solving the linear system 

N i-b 

^t]wj^ x''p{x)dx; n = 0,l,...(Ar- 1). (5) 

j=i 
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Interpolation. The truncated Fourier expansion can be converted into an interpolative formula with the aid of the 
cardinal function Gj{x) with the property Gj{ti) = Sij. This function can be obtained by construction 

G,{x) = 0jv+i(a;)/[0^+i(tj)(a; - t,)] (6) 

where prime denotes the derivative and the interpolative formula reads 

N 

u(x)=^ii(i,)G,(x). (7) 
The cardinal function is a polynomial of order {N — 1) and may be written as a superposition of the basis functions 

N 

G,(x)=^MrV.(x) (8) 

1=1 

where the coefficient matrix can be determined by taking advantage of the orthogonality of the basis. Using 

the Gauss integration formula, we obtain 

Mr^ = 4,,{t,)w,/{c^,,(l)^) (9) 

This matrix can be immediately inverted by applying Gauss integration formula for expressing the orthogonality 
condition for the basis 

N 

<pix) (j)k{x) p{x) dx = 5ik{(t)i,cj)i) = '^<Pi{i]) (10) 

i=i 

and this result is exact. Since the orthogonality condition has here the form M^^ ■ M = 1, it remains to read off the 
the matrix M 

M,j^(^,{t,). (11) 

The function u{x) is determined by supplying the coefficients Cj occurring in the spectral expansion formula but the 
latter is completely equivalent to the interpolative formula that instead of the Cj employs the grid values u{tj). There 
is a linear relation between these two arrays provided by the M matrix which can be written in a concise form 

[u]=M-[c]. (12) 

The array of the coefRcients [c] is obtained by solving the linear system In the semi-spectral method, this gives 

N 

J2{Lcl),{x)}^=t,c,=f{t,), i = l,2,...N (13) 

whereas the Galerkin scheme yields 

N 

Y,{<l>k,Lcl),)c,^{^kJ), k^\,2,...N (14) 

However, if the integration in 114|l is effected by the Gauss method, the resulting equations in both, semi-spectral and 
Galerkin schemes turn out to be identical and the two methods are completely equivalent. Indeed, multiplying both 
sides of (|13|) by (f>k{U) Wi and carrying out a summation over i, we immediately obtain H14|) . In many applications the 
interpolative form ^ might be more convenient, in which case uitj) play the role of the unknowns. The latter are 
determined by solving the linear system 

N 

Y,{LGj{x)}.,=uu{tj) = f{t,), z = l,2,...7V (15) 
i=i 

The resulting u{tj) are subsequently used in (jT)! to obtain the ultimate solution of the problem, as advertised at the 
beginning of this section. 
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In this paper we are going to use the Chebyshev polynomials as the basis functions. There are strong arguments cor- 
roborating this choice. Firstly, the Chebyshev polynomials optimize the interpolative procedure, as will be discussed 
shortly. Secondly, the values of the grid points tj for all N are given by an elementary analytic expression whereas for 
other polynomials they would have to be computed numerically for every N. Thirdly, the Fourier coefficients of the 
function and the appropriate coefficients specifying the derivative are connected by a simple linear transformation. 
Fourthly, the same is true for the antiderivatives. Simplifying matters, one may say that by differentiating or integrat- 
ing a set of Chebyshev polynomials we always obtain a set of Chebyshev polynomials and therefore these operations 
can be reduced to simple matrix manipulations. In the next Section we have presented all the necessary tools for 
solving some typical differential and integral equations occurring in quantum mechanics. Subsequently, we consider 
the application of the semi-spectral Chebyshev method. Due to space limitations some selection was inevitable and 
we decided to confine our interest to the time independent non-relativistic two-body problems. We concentrate our 
attention on the situation when the forces take the form of a short ranged potential considering in sequence the 
continuous and the discrete spectrum. The solutions are obtained both in configuration and in momentum space 
using the appropriate Schrodinger, or the Lippmann-Schwinger (from now on referred to as L-S) equation, with and 
without the Coulomb interaction. These kind of problems are often encountered in nuclear and particle physics and 
we offer a variety of complete solutions. 

In the last part we present some numerically solved examples. The short-range interaction is modelled by three 
different spherically symmetric potentials. The first of them^the exponential potential can be viewed as the simplest 
prototype of the forces rapidly falling off with the distance Q . More complicated shapes can be obtained by taking 
suitable superpositions of exponential potentials with different ranges. Our second choice, the Morse potential 0, is a 
particular combination of two exponential potentials of different ranges where one of them is repulsive while the other 
is attractive. This potential has been very popular Q as a model for vibrational states of diatomic molecules but it has 
been also used in nuclear physics to simulate a repulsive core nucleon-nucleon interaction Q . As the other extreme one 
may take an infinite superposition of exponential potentials forming a geometric series. The explicit summation of this 
series leads to our third choice, the familiar Hulthen potential which close to the origin develops a 1/r singularity, 
like the Coulomb potential, and yet exhibits an exponential fall-off at large distances. The meaning of the Hulthen 
potential goes far beyond academic interest as it has been widely used in nuclear and particle phenomenology 10] llj, 
atomic physics solid state |0, quantum chemistry etc. Here, for us the important fact is that for i = all the 
three potentials mentioned above admit analytic solutions and the situation has been summarized in the Appendix 
[Lj JiLj- In most instances the semi-spectral method turns out to be so accurate that the standard procedure, where 
the convergence is examined by comparison with the results obtained by "other methods", almost certainly would be 
meaningless. Therefore, for calculating the relative errors we always utilized the exact results. 



II. CHEBYSHEV INTERPOLATION 

A. Chebyshev cardinal function 

Let us consider an independent variable t defined in the [—1,1] interval (our conventions and notation in this section 
are consistent with ) The Chebyshev polynomial (t) of the first kind and of degree N is defined 0] by the 
formula 

TN{t) = cos(7Varccos(i)). (16) 
It is evident that Tpf{t) has N zeros in the [—1,1] interval and they are located at the points 

tfc = cos[7r(fc- i)/7V]; k = l,2,...,N (17) 
The Chebyshev polynomials can be also generated from the recurrence relation 

T„(0 = 2ir„_i(t)-r„_2(t) 

with TQ{t) = 1 and Ti{t) — t. The Chebyshev polynomials of the second kind, denoted as Un{t), are obtained from 
the same recurrence relation but with different starting values: U-i{t) — and Uo{t) — 1. 

Usually, two choices of grids are considered: (i) the classical Chebyshev mesh-points given in (|17() , or, (ii) the zeros 
of the derivatives of Tjv(i) which is the same as the zeros of UN-i{t). The choice (ii) is called Lobatto mesh, for it was 
shown by Lobatto that the zeros of the derivative T'^{t) supplemented by the end-points ±1 have all the advantageous 
features of the Chebyshev mesh. Furthermore, if a function is approximated by Chebyshev polynomials supplying the 
boundary values of the function immediately pins down two coefficients. The other attractive feature shows up when 
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the Lobatto mesh is used in quadratures and when N and subsequently 2A'' points are taken to obtain useful and 
computable error estimate. In the bigger mesh half of the points coincide with those of the smaller mesh and therefore 
the function values calculated at the smaller mesh can be later reused reducing thereby the total number of function 
evaluation. The Lobatto mesh, however, is not always preferable and especially for functions which exhibit end-point 
singularities more convenient might be the classical Chebyshev mesh. Apart from personal prejudices, there is little 
difference between these two schemes in terms of accuracy, or convergence • 

The Chebyshev polynomials are orthogonal in the interval [—1, 1] over a weighting factor — and from the 

continuous orthogonality relation a discrete orthogonality relation can be derived. In this work we shall use the 
classical Chebyshev mesh ()17|l and the appropriate orthogonality relations, take the form 



N 



fe=l 



J2 ^^i^k) T,{tk) = 5,, ^(1 + 5^i) (18) 



where i < N and j < N. 

A continuous and bounded function f(t) can be approximated in [—1, 1] interval by the expression 

N 

(19) 

where the primed sigma denotes hereafter a summation in which the first term should be halved. The spectral 
coefficients Cj can be obtained by multiplying H19|l by Ti_i(t) and then setting t equal to tk and, finally, performing 
a summation over k. Making use of l(TH)l . we get 

2 ^ 

^^■ = ]^E^^^fc/(^'«)' (20) 

k=l 

with A/jTj = Tj_i(tfc) where 

T,_i(tfe) = cos[7r(fc - - 1)/N]. (21) 

We may regard l(T^ as a generalized truncated Fourier expansion. Recalling that Tiv(t) vanishes at all t ^ tj, itTD)) 
can be rewritten as 

where prime stands for the derivative. This is nothing else but an interpolating formula, which may be cast to the 
familiar Lagrange form 

N 

i=i 

where we have introduced the Chebyshev cardinal functions Gj{t), with the property Gj{ti) = Sij. Indeed, inserting 
(|2nj) into (O, we obtain 

2 ^ 

^^W = l^E'^-i(^J-)^-i(*)- (24) 



1=1 



It is evident that there are two options for reconstructing the function f(t): either from the expansion (|19|) knowing 
the coefficients Ci, or by supplying the values of f{tj) and using H23|l . These options are completely equivalent and 
knowing the mesh-point values f{tj) of the function, the corresponding Fourier coefficients can be obtained from H2U|I . 

Similar procedure may be applied for a function of two, or more, variables. Thus, a simple extension of H23() to the 
case when the function depends upon two variables x and y, is 



N 



f{x,y)= G,{x) f{U,t,)G,{y). (25) 

Owing to their rapid, often exponential, convergence rate when the number of mesh-points N increases, the Cheby- 
shev interpolation methods have proved to be extremely efficient in solving various problems of mathematical physics. 
This success has a very solid background resting on two celebrated theorems which we quote here without proof 0. 
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Theorem 1 (Cauchy interpolation error theorem). Let f(t) be a function sufficiently smooth so that it has at 
least N + I continuous derivatives on the interval [—1, 1] and let PN{t) be its Lagrangian interpolant of degree N for 
t G [—1, 1]. Then, for any grid ti G [—1, 1] and any t G [—1, 1] the upper bound of the interpolation error is given by 

1 ^ 

f{t) - P^{t) < ^^^^ /(^+^nO - (26) 

for some ^{t) G [a,b]. 

It is wortli noting that the product term in (|26|l is itself a monic polynomial (leading term coefRcient equal unity). 
The above Cauchy theorem indicates that without specifying the function f{t), the only way to minimize the inter- 
polation error is to minimize the product term in (|26|l and this in turn can be achieved by a judicious choice of the 
interpolation points. The answer how to accomplish that brings the following Chebyshev theorem 2]. 

Theorem 2 (Chebyshev minimal amplitude theorem). Out of all monic polynomials of degree N, the unique 
polynomial which has the smallest maximum on [—1,1] is the Chebyshev polynomial 7Ar(t) divided by 2^~^, i.e. all 
monic polynomials Q n {t) of degree N satisfy the inequality 

max\QN{t)\ > max|TAr(i)/2^-i| = 1/2^-^ (27) 

for all t G [-1, 1]. 



B. Anti-derivative 



The optimized interpolation is highly satisfactory but there is much more to come. As we shall see in a moment, 
by using the spectral representation simple operations like differentiation or integration can be reduced to matrix 
multiplication. We start with integration, introducing the anti-derivative of f{t), defined as 

F-(t)= £^ f{x)dx. (28) 
The anti-derivative might be sought in the form of an appropriate Chebyshev expansion 

N 

i^-(i) = ^'QT,_i(t). (29) 

in which case our objective would be to calculate the coefficients Cj occurring in H29|l . If the function in the integrand 
of (|28|) is represented by its Chebyshev expansion H19|l . the integration in (|28|l can be done analytically (Curtis- 
Clenshaw integration [13 )■ This solves the problem because the coefficients Cj in H29() must be expressible in terms 
of Cj but these relations are not particularly simple. El-gendi observed |19| that much simpler relations would be 
obtained if instead of the spectral coefficients the function values at the mesh-points were used |0| [23| • In the latter 
case the array of the mesh values of the anti-derivative F^{ti) would be connected with the array of the integrand 
mesh values f{tj) by a simple linear transformation which bears a universal character being independent upon the 
shape of f{t). Indeed, by inserting (|23|) in (|28|) . we end up with a remarkably simple quadrature rule 

N 

where W^j~ may be viewed as a generalized weighting factor. For a fixed iV, this is a constant matrix, viz. 

Wr^£^G,{t)dt (31) 



which will be obtained in an analytic form. Inserting (|24|) in ()31|l. we get 



^7 = §E' f T,^iit)<^^tn-i{t,) 

k=l -^-1 



(32) 
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and, as seen from (|32|l . the matrix W may be written as a product W = S ■ d ^ ■ M, where M is given in H21|l . 
d = diag(7V, ^N, ^N, ^N), and S", is 

Sr^= jy,-i{t)dt. (33) 

The integral in (|33|l can be obtained from an analytic expression 

[Ti(t,) + 1; forj = l 

dt = } [T2{U) - l]/4; for j = 2 (34) 

[r,(i,)/(2j) - r,_2(t.)/2(j - 2) + (-I)Vj-(j - 2); for j > 2 

where the right hand sides of l|34(l are readily computed with the aid of (|21|l . It is worth noting that the discrete 
orthogonality relation (|18|) can be written in matrix form as M ■ M — d with tilde denoting the transposed matrix. 
Since the matrix M has orthogonal columns, we have M^^ — M d^^ . When the array F~{ti) has been determined, 
the expansion coefficients of the anti-derivative may computed from 

2 ^ 
i=i 

and the anti-derivative F^(t) could be reconstructed at any point either from H29f) . or from ()23|l after making the 
substitution / — > F~ . 

Alternatively, we could have defined the anti-derivative as 

i^+(i)--^ f{x)dx (36) 
in which case the appropriate integration rule, takes the form 

N 

F+{U) = ~J2 ^^/fe) (37) 

where = ■ d~^ ■ M and for S~^, we obtain 

^1 (-T,iU) + l; forj^l 

/ T,^i{t)dt= }-[T,{t,)-l]/A; for J = 2 (38) 

[-TjiU)/i2j) + T,^2{U)/2{j-2)-l/j{j^2); for j > 2 

In the following we shall need both forms of the anti-derivatives. 



C. Integration 

A definite integral in the [— 1 , 1] limits might be calculated in a similar way as the anti-derivative but we wish to 
consider the numerical evaluation of a slightly more general integral 

/[/]= J'^ k{t)f{t)dt (39) 

where k(t) is a real absolutely integrable function, which need not be continuous or of one sign, while f{t) is any 
continuous function. The integration will be performed by using the product integration form 

N 
i=l 

where the weights Wi are determined by requiring the above rule to be exact when f(t) is any polynomial of degree 
< N. This can be regarded as a variation of the Clenshaw-Curtis integration and indeed reduces to that method when 
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the mesh-points are taken to be the Lobatto points an k{t) is set to unity. It has been shown [211 ii'-Sit the integration 
rule H4U|I for iV ^ oo converges to the exact result for all continuous functions f{t), provided k{t) satisfies a rather 
mild integrability condition, namely 

1 

\k{t)\Pdt<oo (41) 

-1 

for some p > 1. Two choices of k{t) are of particular interest. The first is when k{t) is taken to be the Chebyshev 
weighting function k{t) = l/Vl — and we end up with the standard Gauss-Chebyshev quadrature 

1 N 

jit)^f=^-^T.m^ (42) 



with Wj = Ti /N for all j . The second case of practical interest is that in which k{t) = 1 and we have 

/ /(i)dt-5]^,/(i,) 



(43) 



with the Chebyshev weights 

N 



2 

N ■ 



JV .1 

Y,'Tk-i{tj) / Tk-i{t)dt. 
The integral occurring in 1441) can be easily evaluated, viz. 

fju-i{t)dt^[{-lf-l]/[k{k~2)] 
and inserting H45|l in (|44|l . we arrive at the ultimate expression for the weights 



(44) 



(45) 



. {(N-l)/2\ 

= V' (46) 

J TV ^ 4^2 - 1 ^ ' 

i=0 



It can be proved that the weights are all positive and setting f{t) = 1 in we obtain the sum rule 

N 



N 
i=l 



D. Singular integrals 

The Chebyshev quadrature can be extended to comprise also singular integrals. To this end we consider product 
integration (|39|l in which k(t) may contain isolated singularities [22|. Perhaps the most important of them is the 
Cauchy integral, and we wish to establish the quadrature rule for the principal value integral of the form 

l^dt^Y.^,{z)f{t,) (47) 

where the weights ujj{z) depend upon the location of the singularity. It is assumed in the following that |z| < 1 and 
that the function f{t) is free from singularities on the < —1,1 > interval of the real axis. Our goal is to calculate the 
generalized weights ujj(z) and inserting H23|l in the left hand side of (|47|) . we have 

2 ^ 

^J-W = ]vE'^'^-ife)^''-i(^) (48) 
fe=i 
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with /fe_i(z) denoting the principal value integral 

Ik-i{z) 



t - z 



At 



which will be calculated in an analytic form. Indeed, a simple subtraction, yields 



^ Tk-i{t)-Tk-i{z) 1-z 

dt + Tk-i[z) log — — 

1 t — z 1 + z 



(49) 



(50) 



and the integrand on the right hand side of (|50|) being now regular, can be represented as a superposition of Chebyshev 
polynomials 



Tk-i{t) - Tk-i{z) 



fe-2 



Z 



= Uk-2-^{z)T^{t) 



i=0 



rendering the t-integration straightforward. With the aid of (|45|l , we obtain the ultimate expression 

Ik-i{z) = Sk-i{z) + Tk-i{z) log[(l-z)/(l + z)]; 

with 

Sk-i{z) = -ZV"' ^ 1^ ~\ ^ Uk-2-i{z) 



(51) 



(52) 



(53) 



and H52|) - (|53|) used in H48|) allows to calculate explicitly the weights. The presence of the end-point logarithmic 
singularities is apparent in H52|l and at these values the principal value integral (|49() becomes undefined. It is interesting 
to note that the relation between the function /„(z) and T„(z) is analogous to that between the Legendre function 
Qn{z) and the Legendre polynomial Pn{z). It will be convenient to separate the singular term writing the weights in 
the form 



u;,iz) = Co,{z) + G,{z) log [(1 - z)/(l + z)]; 
where Gj{z) is the cardinal function and the regular part, given by 

2 ^ 



(54) 



(55) 



k=l 



is just a polynomial in z. The ultimate quadrature rule for a general integral of the Cauchy type, reads 



' fit) dt 

_i t — z^ie 



— z 

j=i J 

J2^,{z)f{tj)+f{z) hog^iiTT 
j = l 



\z\>l 



(56) 



\z\<l 



Setting f{t) = l in (|47|l . a sum rule is obtained 



N 



^Wj(z) ^ log 



1 - z 
1 + z 



providing a convenient check of the computed weights. 

The other important case to be considered here is that of the logarithmic singularity. The appropriate automatic 
quadrature, takes the form 



J J{t) log\t-z\dt^< 



( N 

^u;j/(ij)log|tj -z|, |z|>l 

N 

Y,^,{z)f{t,), \z\<l 



(57) 
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and our goal is to determine the weights ^j{z). The assumptions concerning f{t) and z remain the same as in the 
Cauchy integral case (|47|) . In particular, for the time being, we assume that |z| 1 but eventually this restriction will 
be lifted. Upon inserting H23|l in the left hand side of 1)57(1 . the dedicated weights ^j{z) can be written in the form 

2 ^ 

^A^)^J^T.'T^-^(^^)J''~^(') (58) 
fe=i 

with Jfc_i(z) denoting the integral 

Jk-i{z)^J Tk-i{t) log\t- z\dt, (59) 

which can be evaluated analytically. In order to do that, we take advantage of the identity which holds for 
non-negative n 

{\n{z), n = l. 

with prime denoting the derivative. Inserting the above identity in H59I) . the integration by parts, gives 

l0g|l-z|-(-l)'=l0g|l+z| Ik[z) , I\k-2\{z) 



Jfc_l(z) = 



k{k-2) 2k 2{k~2' 

1 + z 



., k^2 



(60) 



where k — 1,2, - • • ,N. Substituting the explicit forms of Jk-i{z) in ((58() completes the derivation of the weighting 
functions flj{z). 

The integral containing a logarithmic singularity H57|l . unlike the Cauchy integral, does exist even when the singu- 
larity coincides with either of the integration end points. Therefore, all Jk-i{z) in H60|l go, as they should, to a finite 
limit when z ^ ±1 and the explicit expressions, are 

'21og2-2, fc = l 

kijr^) M + 2or:iy 

where 5'fc(±l) and S'fc_2(±l) are computed from H53I) . With (|61|l in hand, we can evaluate rij(±l) from 158|1 . 
establishing the automatic quadrature (|57|l for the case when the singularity is located at the integration end points. 
Setting f{t) = l in ((57|l . we arrive at the sum rule 

N 

^A^) - Mz) = (1 - 2) log (1 - z) + (1 + z) log (1 + z) - 2 

j=i 

which might be useful for checking purposes. The automatic quadrature just presented has usually a very high rate 
of convergence and compares favorably with a direct method of calculating the integral (|57() where the singularity is 
eliminated by subtraction. 



E. Differentiation 



The Chebyshev expansion is also well suited for performing differentiation. In fact, we could invert our formulae 
for integration to derive the appropriate expressions for differentiation. It is better, however, to provide an ab initio 
derivation. Upon differentiating (|19|l . we have 

am _^ dr,_i(t) 
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but this function has also a standard Chebyshev expansion 

N 



dfit) _ „ 



/'(t)=^'c;.T,_i(0 (63) 



At 

where denote the appropriate Fourier-Chebyshev expansion coefficients. Clearly, denote the appropriate are 
related to the expansion coefficients for the function Cj. The same is true for the f'{ti) and f{ti) and there is a 
linear transformation connecting these two arrays. To find the transformation matrix explicitly we set in (|62(l t = ti 
inserting for Cj formula (|2(J|) . As a result, we obtain 

N 

where the matrix Dij is 

2 ^ 

D^ = j:^Y.'Tk-i{U)n-i{tj). (65) 

fc=l 

The derivative of the Chebyshev polynomial occurring in l|65ll can be can be easily computed as 

TU{U)^{k-l)Uu-2{h) (66) 

where Uk~2{tj) denotes a Chebyshev polynomial of the second kind. 

When the array f {ti) is known, the expansion coefficients are computed from 

2 ^ 

'^^^nH ^^^^ ^'(*^) (67) 



N 



and the fmiction f'{t) can obtained form (|63|) . 



F. Arbitrary interval [a, b] 

In the above considerations the independent variable t was confined to the [—1,1] interval. This restriction may be 
lifted by introducing a new independent variable x defined in an arbitrary interval [a, h]. A linear mapping generates 
the appropriate Chebyshev mesh in [a, h\ 

Xj=\{h + a)^\{h-a)tj] j = 1,2,3, ...,iV (68) 

The cardinal functions (|24|l constituting the basis for the interpolation formula depend now upon the argument 
obtained by the inverse transformation. The extension of H23|l . is 

fix) = /(-,) (69) 

j=i ^ / 

and the Gauss-Chebyshev quadrature reads 

,fc N 



/ /(x)da; = i(6-a)^Wj/(a 

Jo. ,-=1 



(70) 



with the weights given in H46|l. The modifications for other automatic integrations are rather obvious. 

The anti-derivative is obtained in a similar way: the weights acquire the ^{b — a) factor, induced by the change of 
variable, and the integrand must be evaluated at the mesh-points Thus, the extension of is 

F-{x,)= f{x)dx = ^{b-a)Ywrf{x,). (71) 
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where the weighting matrix W~ once and for all is given in (|32() . 

In order to apply the Chebyshev method for calculating semi-definite integrals there are two possibilities: (i) trun- 
cation followed by a linear mapping, or (ii) non-linear mapping. Truncation method replaces the infinite integration 
domain [0, oo] by a large albeit finite domain [0, R\. Non-linear mapping consists in a change of variables devised in 
such a way that the original integral is converted into an integral which is in the [—1, 1] limits. Clearly, there would 
be a countless number of such transformations but the simplest one is the rational mapping 

r = R{l + t)/{l-t), (72) 

which relates the original variable r G [0,oo] with a new variable t G [—1, 1]. Dealing with improper integrals has 
its price: in both cases an adjustable parameter R appears, in one case as a cut-off in the other (|72l) as a scaling, 
or a slope parameter, whose value could be optimized at the expense of some trial-and-error procedure. In general, 
the truncation method is more sensitive to the variation of R. Actually, the two numerical parameters N and R are 
correlated and some experimentation is inevitable if one wants to minimize the error by increasing them simultaneously. 
In addition to the appearance of the truncation or scaling parameter R the infinite integration range induces also 
another unwelcome feature in that the increase of R brings all the poles and branch point singularities of the integrand 
(or its derivatives) closer to the real axis after the original variable has been rescaled to the [—1,1] interval. The 
presence of singularities close to the integration region in consequence leads to a deterioration of the convergence rate. 



III. APPLICATIONS 



The implementation of the semi-spectral Chebyshev method is very simple because all that is needed is a standard 
linear system solver plus a set of procedures calculating the cardinal functions Gj{t), the anti-derivative matrices 
W^, the differentiation matrix D and the dedicated weights Wj ,u!j{z),il,j{z). Explicit means for getting all of them 
have been described in detail in Sec. II. Assuming then that the toolbox specified above is available, we are going 
to review some applications of the semi-spectral method that are of relevance to quantum mechanics. The scheme 
remains the same in all problems: the equation to be solved is first discretized and on the grid the integration and 
differentiation operators are represented by appropriate matrices. As a result, the mesh-point values of the sought for 
function can be pined down by solving a linear system of algebraic equations. Finally, by feeding the interpolation 
formula H23() with the mesh-point function values, the solution at any point can be obtained. 



A. Second order differential equation 



The second order differential equation of considerable interest in physics (Schrodinger equation), has the form 

y"{x) + p{x) y{x) — q{x); a < x <b (73) 

where p{x) and q{x) are given functions. The unknown function y{x) and its derivative are assumed to satisfy some 
boundary conditions at a or/and at b depending on the nature of the physical problem at issue. We consider first the 
case when the values of y{a) and y'(a) have been assigned. When (|73|l is integrated twice from a to a;, we obtain 

y{x)-y{a)-y'{a){x-a)+ As dtp{t)y{t)^ ds dtq{t). (74) 



Introducing the Chebyshev mesh (|58|l for the external variable x and performing Gauss-Chebyshev integration over 
the internal variables s and t, (|74|) takes the form of a matrix equation 

{l + l{b-a)^W- ■W-o[p]}.[y]^ 
^y{a)[l]+y'{a)[x-a] + \{b~arW- -W- -[q] 

We use the notation where a vector [/] contains the values of the function f{x) at the mesh points (|68|) . i.e. [/] 
stands for the array [f(xi), f{x2), ■ ■ ■ , /(xjv)]^- The symbol and [A o B]ij = AijBij denotes the Schur product of the 
matrices A and B and 1 is a unit matrix. Writing the Schrodinger equation (|73|) in the form H75(l . the problem has 
been reduced to that of solving a linear system easily handled by standard methods. When the vector [y] has been 
obtained by solving the system (|75|l . the function y{x) can be reconstructed at any point x S [a, b] with the aid of the 
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interpolative formula H69|l . Knowing y{x), the derivative y'{x) can be calculated either by Chebyshev differentiation 
(|64|l or from the integral expression 

y'{x)=y'{a)~ r Atp{t)v{t)+ C dtq{t) (76) 



which after discretization involves only a matrix times vector multiplication. The derivative is then obtained by 
interpolation. It has to be noted that in the above considerations the interval [a, 6] has been arbitrary. This means 
that eq. (|73|l may be solved using composite Chebyshev integration that introduces a partition of the interval 
[a, 6]. Performing integration in the first subinterval provides us with the values of y{x) and y'{x) at the end-point. 
Subsequently, they can be used to specify the initial conditions for the integration in the second subinterval, and so 
on until integration in the last partition has been accomplished. This option might be particularly useful when the 
functions p{x) and q{x) are rapidly varying and their representation would require exceedingly large Chebyshev mesh. 
Composite integration is in such a case the simplest remedy. 

With different boundary conditions the procedure is similar. Often we know y{a) and y' (b) in which case eq. H73|l 
is first integrated from x io h and next from a to x. As a result, we obtain 

/■a: ph px ph 

y{x)-y{a)-y'{h)[x-a)~ As dtp{t)y{t) = - ds dtq{t). (77) 

J a J s J a J s 

Putting the external variable x on the Chebyshev mesh and using Gauss-Chebyshev quadratures, l(77|l takes the 
form 

{l-'^ib-afW- ■W+o[p]}.[y] = 

= y{a)[l]+y'{b)[x^a]~\{b^arW--W+-[q]. 

It is evident that also in this case we end up with a linear system which, naturally, is different from (|75l) . In either 
case, when the vector [y] has been determined the ultimate solution of (|73|l is given by the interpolative formula l|69|l . 

B. Volterra integral equation of the second Itind 

We write the Volterra integral equation of the second kind, as 

y{x)-X f k{x,s)y{s)ds^f{x) (79) 



where a < s,x < b and A is a given parameter. Both, the function f{x) and the kernel k{x, s) are also assumed to 
be provided, whereas the function y{x) is the sought for solution of H79|) . Using Gauss-Chebyshev quadrature in H79|l . 
we are led to a system of N linear algebraic equations 

{l-X^{b^a)KoW-}-[y] = [f]. (80) 

The solution [y] of the system H80II used in the interpolative formula (|69|l yields y{x) at an arbitrary point from the 
[a, b] interval. 

When the independent variable x appears instead as the lower integration limit, i.e. if H79|) is replaced by the 
Volterra equation 

yix)-X [ k{x,s)y{s)ds = f{x) (81) 

J X 

the procedure is similar and the only change required is W~ — > W"*". In the end we arrive at the system of equations 

{l-\\{b~a)KoW+] ■[y]^[f]. (82) 
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C. Fredholm integral equation of the second kind 

The Fredhohii integral equation of the second kind can be written, as 

y(x)-A / k{x,s)yis)ds^ f{x) (83) 

J a 

where f{x) and k{x, s) are known functions assumed to be regular in [a, b] whereas y{x) is the unknown function. The 
Chebyshev method is formally akin to the Nystrom method in which Gauss-Chebyshev quadrature has been adopted. 
Similarly as in the Volterra eq. case just considered, we end up with a system of linear equations 

{l^X^ib-a)Kow}[y] = [f]. (84) 

where [K o w]ij = k{xi, Xj) Wj with the weights Wj given by H46|) . 

In many physical applications the kernel fc(a;, s) in (|83|l might be continuous on the diagonal x — s but the derivative 
of the kernel exhibits a discontinuity at this point. A typical situation occurs when 

\ki{x,s) tors<x 

HX,S) = S , . ^ f . (85) 

k2{x, s) tor s > X 

and if fci and fc2 are different functions, the kernel is referred to as semi-continuous. The Chebyshev method is well 
suited to handle such situation and the appropriate extension of (|84|l . reads 

{1 - A i(6 - a) o W + ^(2) o w+) }-[y] = [/] (86) 

where K^^^ — ki{xi,Xj) and = k2{xi,x.j). 



D. Singular integral equation 

The singular integral equation with a Cauchy type singularity, can be written as 

y(a;)_A rM^y(,)ds = /(a;) (87) 
J a s~ z 

where z is a real parameter and the integral is defined in the principal value sense. When z falls outside the integration 
limits the integral is not singular and eq. H87I) becomes a Fredholm equation considered above. Therefore, we confine 
our attention to the case when z lies within the integration limits {b > z > a). Apart from the pole at s = z the 
kernel is assumed to be regular on the real axis in the integration range. The singular quadrature rule introduced in 
(|47|l allows to apply the Chebyshev method in the same way as in the Fredholm equation case. This procedure may 
be viewed as an extension of the Nystrom approach and we end up with a system of linear algebraic equations 

{l-Al^oc.(r)}.[y] = [/]. (88) 

where \K o w(r)]ij = k{xi,Xj) ujjir) and t = [2z — {b + a)]/{b— a). Formally, the solution of eq. H87|) and the solution 
of (|83|l look very much the same and what differs them are the weighting factors: the singular quadrature employs 
dedicated z-dependent weights Wj(T) given in H48|l . For a singular integral equation, however, the Fredholm alternative 
does not hold and the particular solution of the inhomogeneous equation obtained above has to be supplemented by 
the general solution of the homogeneous equation which must be obtained by a different method. 

The weakly singular integral equations can be handled in a similar manner. We shall mention briefly the case when 
the kernel has a logarithmic singularity 

y{x)^\f k{x,s)log\s-z\yis)ds^fix) (89) 

J a 

where a < z < b and the function k{x, s) is taken to be regular within the integration range. The integral in (|89|1 can 
be approximated by the automatic quadrature (|57l) and we are led to a linear system of algebraic equations 

{ 1 - Ai (6 - a) K o [u; log i (6 - a) + r! (r)] } • [y] = [/] . (90) 

where {K o [w log ^{b— a) + i}{T)]}ij = k{xi,Xj)[wj log i (6 — a) + r2j(r)]. Here are the Gauss-Chebyshev weights 
and the dedicated weights ^^^(t) are given in lt55|) with r = [2z — {b + a)]/ {b — a). 

An example of H87|) is the Omnes equation 23] which arises in the dispersion theory of final state interaction. The 
equation H89|) appears in the Faddeev approach to a three-body problem Although the above results could be 
immediately applied in these problems, a discussion of these topics is beyond the scope of the present paper. 
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E. Integro-differential equation 



The Chebyshev method is also very well suited to handle integro-differential equations. As an example we may 
consider the equation 

y' (x) + p{x) y{x) — q{x) + / k{x, s) y{s) ds; a<x<b (91) 

J a 

where p{x),q{x) and k{x,s) are given functions. We also assume that the value of y{a) has been assigned. Upon 
integration of H91|l in the limits from a to x, we obtain 

y{x)~y{a)+ p{s)y{s)ds^ g(s) ds + / / k{t,s)y{s)dtds. (92) 

•/a -/a J a -J a 

Introducing Chebyshev mesh, we are led to a Hnear system 

{1 + i(6 - a)W- o [p] - i(6 - a)W- ■Ko[w]]-[y]^ y{a)[l] + \{b - a)W- ■ [q]. (93) 

Another type of integro-differential equation that might be encountered in physical applications is the Schrodniger 
equation with exchange interaction. Consider then a second order homogeneous equation 

y" (x) + p{x) y{x) = I k{x,s)y{s)ds] a<x<b (94) 



where p{x), q{x) and k{x, s) are given functions. We also assume that y{a) and y'{a) take assigned values. When 
is integrated twice from a to a;, we obtain 

y{x) - y{a) - y' {a){x - a) + ds dtp{t)y{t)^ / k{t,s)y{s)dtds. (95) 



Upon introducing the Chebyshev mesh, we end up with a liner system of equations 

{1 _ 1 (5 _ af W W o [p] ~\{b-afW--Ko[w]]- [y] = 

^y{a)[l]+y'{a)[x^a]. ^ ' 

In the Chebyshev approach the inclusion of exchange forces brings only a minor complication. The amount of labor 
required to obtain the solution of the integro-differential equation H94|l is about the same as in the case of the differential 
equation l|73|l . 



IV. QUANTUM MECHANICAL TWO-BODY PROBLEM 



A. Configuration space 



We wish to examine the effectiveness of the pseudo-spectral Chebyshev method by solving explicitly a number of 
scattering and bound states problems but before we do that, we need to establish our notation and assemble some 
well known theoretical tools. For the sake of clarity, we confine our attention to a single channel, time independent 
two-body problem. To simplify matters even further, we assume that the "particles" are deprived of any internal 
degrees of freedom. The adopted interaction has a form of a local, short-ranged, spherically symmetric potential V{r) 
which depends only upon the mutual separation r between the particles. The shape of V{r) is arbitrary and a singular 
behavior when r — > is not excluded with the restriction that this singularity would be not worse than 1/r. Since 
the 1/r singularity requires special care, at the top of V{r) we introduce explicitly a point charge Coulomb potential 
Vc(r) — aZ/r where a is the fine structure constant and Z is the charge number. 

The radial part of the appropriate Schrodinger equation, is 

dV(p,?')/dr2 + [p^ -£{e+l)/r^ - 2^Vb(r)] ■ipip,r) = 2fiV{r) ■iP{p,r), (97) 

where p, fj, and £ denote, respectively, the center-of-mass momentum, the reduced mass and the orbital momentum. 
Since in (|97|l different £ do not mix, the label £ on the wave function il){p,r) has been dropped. The regular solution 
of H97|) is defined by the boundary condition 

i/'(p,r)-/+i for r^O. (98) 
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In a scattering problem the momentum p takes real non-negative values whereas in a bound state problem the 
momentum will be imaginary. We shall deal with both problems in sequence commencing with the continuous 
spectrum case where our objective is to calculate the scattering phase shift. 

In absence of V{r), the two linearly independent solutions of l|97|) . and referred to in the following as the "free 
solutions", are 

f{p) = Fe{v,p), 
g{p) = Gi{r],p), 

where Fi{rj^p) and Gi{'i],p) are the standard Coulomb wave functions defined in 4], p ~ pr and rj ~ fj,aZ/p is the 
familiar Sommerfeld parameter. When both, the Coulomb interaction is switched off {Z — 0), and V{r) — 0, the free 
solutions simplify to the form 

f{p) = pii{p), 

9{p) = -pntip), 

where ji and ng denote the usual spherical Bessel functions. 

For r > R where R is some cutoff radius which is sufficiently large as compared with the range of the potential 
V{r), the potential term on the right hand side in H97|l gives negligible contribution and 'ip{p,r) is a superposition of 
the free solutions 

ip{p, r) « A[f{pr) + ianS g{pr)], for r > R, (99) 

where A is a constant amplitude and 5 denotes the phase shift (suppressing the £ label). Given the logarithmic 
derivative of il^ij), R), the phase shift is computed from 

,^ f{pR)^P'ip,R)/^P{p,R)-npR) 

g[pR) ^'{p, R)/i;{p, R) - g'{pR) ^ 

where prime denotes the derivative with respect to R. The free solutions can be also utilized in constructing the 
standing-wave Green's function Gp{r,r') 

Gp(r,r') = -il/p)f{pr^)g{pr,) (101) 

where r< = min(r, r') and r^. — max(r, r'). Because the Green's function H1U1|) is a solution of the inhomogeneous 
wave equation 

d^Gp{r, r')/dr2 + [p^ - £{£ + l)/r^ - 2^ Vcir)] Gp{r, r') = 5{r - r'), (102) 
the Schrodinger equation H97() may be immediately converted into the L-S equation 

/•oo 

^{p,r)^ f{pr) + 2p / Gp{r,r')V{r')^{p,r')dr'. (103) 

"'0 

This integral equation incorporates the boundary conditions (|98|l and (|99|l . Since the incident wave has amplitude 
equal to one, given the solution of p()3|l . the phase shift is obtained from the formula 

tan(5 = -— / }{pr)V{r)il){p,r)dr. (104) 



P Jo 

It will be advantageous casting the equation into a Volterra equation. Indeed, the Fredholm equation H103() can be 
rewritten, as 

V'(p, r) = Cfipr) r [f{pr') g{pr) - g(pr') f{pr)] V{r') V'(p, r') dr' (105) 

P Jo 

where 

C=l-^/ g{pr)V{r)'iP{p,r)dr (106) 
P Jo 
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is a constant which does not depend on r. At first sight, (|105|l involving an unknown constant C does not appear to 
be of much use, but in fact it is possible to get rid of C and eventually pin down its value. Indeed, introducing a new 
wave function function u{p, r), defined by the formula ip{p, r) — C u{p, r), eq. H105(l reduces to a Volterra equation of 
the second kind for the function u{p, r) and this equation is free from unknown constants 

u{p, r) = f{pr) / [f{pr') g{pr) - g{pr') f{pr)] V{r') u{p, r') dr'. (107) 



P Jo 

Solving eq. H107|l for u{p, r), the value of C can be determined by first substituting = Cu in H10t)|) and then solving 
the resulting equation for C. We get 

C==(l + — / g{pr)V{r)u{p,r)dr\ (108) 



P Jo 

and the ultimate formula for the phase shift obtained by substituting ijj ~ Cu in 1104(1 . is 



ta.nd=-C—l f{pr)V{r)u{p,r)dr, (109) 



P Jo 

with C given in H108|l . This formula makes reference only to the solution of eq. H107() . The possibility that the 
configuration space non-relativistic scattering problem might be formulated in terms of an inhomogeneous Volterra 
equation of the second kind was first noted by Drukarev [23 more than half a century ago. Although his method 
soon received the necessary mathematical background , the Volterra equation approach went into oblivion to be 
revived only recently i^27j 28] 29] . 

The above scheme based on the Volterra integral equation may be also applied for calculating the scattering length 
but the calculation of the Coulomb corrected scattering length is slightly more complicated and for this purpose one 
needs a different set of Coulomb wave functions which are analytic and entire functions of p^. Such functions, denoted 
here as ^eip, f) and Qi{p, r), are particular superpositions of Fg and G^, and are available in the literature. Following 
[sof , for p = the regular function takes the form 

([i2e + iy./f3'+^]^i2i+i{2^) z>o 

$^(0,r) = <^ Z = (110) 
[[(2^+l)!//3^+i]V^ J2^+i(2V^) Z<0 

where (3 = 2^a\Z\, and (J2£+i, ^2£+i) are Bessel and modified Bessel functions, respectively The irregular function 
61(0, r) is given in terms of the Neuman Y21+1 and modified Neuman function -K'2^+1, respectively 

'2[/37(2£+l)!]V^X2£+i(2V^) Z>0 
9^(0, r) ^ { r-^/{2i +1) Z = (111) 

-^[/37(2^+l)!]V^y2£+i(2V^) Z<Q 

and with the adopted normalization the Wronskian M^[8£, ^i] is equal unity. The Coulomb corrected scattering length 
Ac is then obtained |30j] from the formula 



Ac = -2t^J $o(0,r)y(r)(/.(r)dr |l + 2Aiy 60(0, r) V (r) 0(r) drj , (112) 
where the function (j){r) is a solution of the Volterra equation 

^{r) - $0(0, r)-2^l ( [$o(0, r') 60(0, r) - $o(0, r) 60(0, /)] V{r') </.(r') dr'. (113) 



The above scheme remains valid in absence of Coulomb interaction. 

Although the method outlined above does not seem to have been tried for bound state calculations, such extension 
is perfectly feasible. Bound states are identified with poles of the T-matrix that are located on the positive part of the 
imaginary axis in the complex momentum plane. Therefore, formula (|109|l must be replaced by the appropriate T- 
matrix expression in which the physical momentum p will be analytically continued to the complex momentum plane k. 
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Using the T-matrix scheme for scattering requires the outgoing wave Green's function Gp^"* (r, r') = — ft(pr> ) f{pr^ ) /p 
where h(x) — g{x) + if{x). The corresponding wave function 4'~^{p, r) is obtained as the sohition of the equation 

/>oo 

i;+(p,r) = f{kr) + 2n Gl,+\r,r')V{r')^+(p,r')dr'. (114) 



Passing to the asymptotic hmit r ^ oo in (|114|l immediately yields the on-shcU T-matrix formula 

T{p) ^ e^^ smd = -— [ f{pr)V{r)'ilj+{p,r)dr. (115) 



P Jo 

Similarly as before, eq. (jll4|l can be converted to a Volterra equation 

^+ip,r) =Df{pr) - ^ / [f[pr')g[pr)-9[pr')f{pr)] V(r')^+{p,r')Ar' (116) 
P Jo 

which has the same kernel as (|105|) but the constant D is different from C and reads 

D = l-— / h{pr)V{r)'iP+{p,r)dr. (117) 



P Jo 

The unknown constant D may be removed by setting ip~^{p,r) — Du{p,r) where the function u{p,r) is a solution of 
(|107|l . the very same equation which was used in the calculation of tan (5. Eliminating tp~^ in favor of u in H115() . the 
ultimate T-matrix formula takes the form 



r(p) = -^/ f(pr)V(r)u(p,r)dr + — [ h{pr)V(j 
P Jo V P Jo 



r)u{p,r)dr\ . (118) 



The T-matrix (|118|) satisfies elastic unitarity constraint 3T^^ = —1. The expression occurring in the denominator of 
pi8|l is recognized as the Fredholm determinant whose zeros located on the positive imaginary axis in the complex 
k-plane would be identified with bound states. Therefore, we need to make an analytic continuation from the physical 
region onto the imaginary axis p ^ in where k is real and non-negative. As a result, the Coulomb wave functions 
and H'^ will be modified, and we have 

Ft{r),pr) ^ i^+i exp(ii^?7) />r), (119) 
Hii'n.pr) -> i-^ cxp(-ii7rry) /i(Kr), (120) 



where fj = fj,aZ/K. The new functions / and h are both real and are given in analytic form 

finr) = C{2nrY+^ e"'^'' M{i+l + fj,2e + 2,2nr); (121) 
h{Kr) = {2KrY+'^ e-'^'^ U{e+l + ii,2e + 2,2Kr), (122) 

with C = \r{£ -|- 1 -I- fi)\/[2{2£ + 1)1] where M{a,b,x) and U{a,b,x) denote, respectively, the Kummer and Tricomi 
confluent hypergeometric functions Q in which all arguments are real (numerical methods of calculating negative 
energy Coulomb wave functions have been presented in [IJ,!^^ and [s^)- In the charge- less case {Z = 0), the above 
functions simplify to the form 

f{Kr)^inr)ifXKr); (123) 
h{Kr) ^ {2/7r){Kr)keiKr), (124) 

where ii{x) and ki{x) are the modified spherical Bessel functions defined in 4]. The adopted normalization is such 
that the Wronskian h{x) f'{x) — f{x) h'{x) is equal to unity. 

The kernel in H105() analytically continued onto the imaginary axis remains real and we infer that u(iK, r) acquires 
merely a constant phase factor equal i^+^ exp (i ^ttt)). Therefore, we set u(iK, r) — i^"*"^ exp (i ^irff) u{k, r) with real u. 
On the imaginary axis the Fredholm determinant A(k) is real, and we have 

A(k) = 1 + — / h{Kr)V{r)u(K,r)dr, (125) 

K Jn 
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where u is the solution of the Volterra equation involving only real quantities 

u{K,r) = /{nr) / [f{Kr')h{Kr) — h{Kr')f{Kr)]V{r')u{K,r')dr'. (126) 

Jo 

The bound states are located by finding the roots of the equation 

A(k) = 0. (127) 

Denoting such root as kq the corresponding binding energy B will be obtained as B = — Kq/2/j,. 

In the literature, the usual method of locating the binding energies utilizes the homogeneous equation (|114|l . The 
configuration space is rather unwieldy for this purpose because the Green's function exhibits a cusp behavior aX r — r' , 
which in general leads to a loss of accuracy. Therefore, momentum space methods have been given preference and 
the task of solving the homogeneous equation reduces to that of an algebraic eigenvalue problem easily handled by 
standard methods. Although, the pseudo-spectral method can cope with such a situation when the kernel is semi- 
continuous l|86|) . but owing to the presence of the matrices the resulting eigenvalue problem in configuration 
space would involve non-symmetric matrices. Unfortunately, the standard methods of solving eigenvalue problems 
are in this case much less accurate than for symmetric matrices, rendering the configuration space method a less 
attractive option. Thus, in the homogeneous equation case the momentum space approach is still preferable. 



B. Momentum space 



In momentum space the partial wave Schrodinger equation, takes the form 

2 roa 

k^(f>{k) + - / Ue{k,k')(l){k')k'^dk' = 2piE(l){k), 



TT 







where Ui(k, k') denotes the l-th wave projection of the local potential 2iiV{r) 

Ui{k',k)^2fi ji{k'r)V{r)Mkr)r^dr. (128) 
Setting E — — K^/2/i for a bound state problem and introducing a new unknown function 



u{k) = k \Jk'^ k2 0(fc) 
we obtain a homogeneous Fredholm integral equation 

/•OO 

u{k) = / IC{k,k')u{k')dk', 
Jo 

with a manifestly symmetric kernel 



/C(fc, k') = -{2/tt) (k/y/k^ + K^) Ue{k, k') (fc7v/fc'2 + k2). 

Following the customary procedure, this kernel may be multiplied by a constant "eigenvalue" A which for an assigned 
K is in fact a function A(k). At a particular value n — for which A(ko) = 1 we would have a bound state. In 
practice, after discretization we end up with an algebraic eigenvalue problem in which k is a parameter. 

The solution to the non-relativistic potential theory scattering problem is obtained from the L-S equation for the 
t-matrix. It is much easier to deal with real quantities and therefore a convenient starting point is the partial wave 
L-S equation for the off-shell K-matrix 



2 n'^dn 

{k'\K,{p')\k)^Ue{k',k)-- u,{k',q)^^{q\Ke{p')\k) (129) 



which is valid for spherically symmetric potentials. When eq. (|129|l has been solved, the t-matrix can be recovered 
from the off-shell unitarity constraint and the resulting fully off-shell t-matrix reads 

ik'\T,ip^±ie)\k) . ik'\K,ip^)\k)^i/J^^^^i^^ (130) 

l±ip {p\Ki{p^)\p) 
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Finally, the on-shell K- matrix is related to the phase shift Si{p), by 

{p\Ki{p^)\p) = -p-' tanSeip). (131) 

The presence of the Coulomb potential creates the well known conceptual difficulties, and e.g., the Coulomb scattering 
amplitude becomes singular on-shell. For {V + Vc), a two-potential problem the scattering amplitude can be written 
rigorously as a sum of the purely Coulomb amplitude plus the Coulomb corrected amplitude associated with the 
short-ranged potential. To escape the difficulties inherent in the ab initio calculation of the Coulomb amplitude in 
momentum space, the latter is usually regarded as known. When this is accepted, the remaining amplitude can be 
obtained quite easily from an extended L-S equation in which the Coulomb Green's function replaces the free Green's 
function. For a repulsive Coulomb interaction the spectral representation of the Coulomb Green's function contains 
only the continuous spectrum and the L-S equation for the Coulomb corrected T-matrix retains its standard form 

{k'\T,{p')\k)^ = U,{k',k),-- / U,{k\q), 2 2- i<l\TeiP^)\k)c 

where the subscript c indicates that Coulomb representation has been used. In particular, the potential term in this 
representation is given as 

/•oo 

Ui{k\k)c^ {2fi/kk') Fi{7^k',k'r)V{r) Fi(r^k,kr)dr (132) 
Jo 

and involves explicitly the Coulomb wave functions with rjk — fJ,aZ/k and rjk' — iiaZ/k' . The T-matrix on-shell 
yields directly the Coulomb corrected phase shift 

{p\Teip^)\p), - -p-' exp(i5|(p)) sin^Kp)- 



V. SOLUTION OF THE SCHRODINGER EQUATION 



A. Continuous spectrum 

In principle, the scattering problem based on the Schrodinger equation could be solved by directly applying H75(l 
where the appropriate boundary conditions are y{a) — and y'{a) = 1. However, there is a better alternative in 
which the proper threshold behavior of the wave function is enforced from the onset by setting ip{Pj = r^'^^(j>{p, r) 
where the unknown function (f>{p, r) is a solution of the differential equation 

d20(p,r)/dr2 = - 2tiV{r)](j){p,r) - [(£ -|- l)/r] d(/)(p, r)/dr. (133) 

Following the well known procedure, the derivative of the wave function with respect to r would be regarded as a 
second function to be determined putting ^(Pi^) = 4>'{PTf)- The wave equation H133|l . takes then the matrix form 



0'(p,r)\ / 1 \((l^{p.r) 

X'{p, r)J y -[P^- 2fi V{r)] ~2(£ + l)/rj {^(p, r) 



(134) 



The system of equations (|134ll will be solved using the pseudo-spectral method by introducing the Chebyshev mesh 
Xj — 5-R[l + cos(7r(j — ^)/N)] in the interval [0, i?]. In order to obtain a system of linear equations for the [0] and 
[x] vectors, we need to impose the boundary conditions at the origin. In view of the fact that the centrifugal barrier 
factor has been already taken care of, the wave function (j) goes to a constant at the origin. Since in the calculation of 
the phase shift we need only the logarithmic derivative of 0, a common normalizing factor in (j) and 4>' is irrelevant, 
and we set (j){p,0) = 1 and x(PiO) = c where c is a constant. For potentials which are less singular than 1/r, setting 
c = would be sufficient to neutralize the singular behavior of the {£ + l)/?" term. The ultimate choice, comprising 
the case when V{r) exhibits a 1/r singularity, is c = /i \imr-^Q rV (r) / {£ + 1). Integrating p34|l on the Chebyshev 
mesh, the resulting system of 2N algebraic equations can be written in matrix form 



.^RW-o[p] 1 + ^RW- o[q]J \[x 



(135) 



where [p]i = {p^ — 2^V{[x]i)} and [q]i = 2(1 + l)/[x],;. When the system (|135|l has been solved, the wave function 
and its derivative are provided by the interpolative formula 

JV 

0(P,O =E' G,{2r/R~l)^{p,x,), (136) 



21 



and 

N 

Xip,r)^j2' G,i2r/R~-l)x{P:X,), (137) 
respectively. Knowing the logarithmic derivative of (j){p,R), the phase shift is obtained from 

. _ fjp^) [(^ + 1)/-^ + x{p, R)I4>{P, R)] - f'ipR) 

g{pR)[ie+l)/R + x{p,R)/Hp,R)]-9'{pR)' ^ ' 

where prime denotes derivative with respect to r. 

The above scheme may be used to calculate the scattering length by going to the limit p 0. Dividing both sides 
of 1138|l by p and subsequently by letting p go to zero, we obtain the scattering length {£ — 0) 

A= -RC{R)/[1 + C{R)], (139) 

where C{R) — i?x(Oj R)/4>{^j R) ^^'^ </'(0i R) ^^^d x(Oi R) ^re the solutions of the system (|134|) in which the momentum 
p has been set equal to zero. 

When the potential V{r) is given as a sum of a short ranged "nuclear" potential plus a point charge Coulomb 
potential, formula 1)1 38(1 remains valid and the resulting phase shift becomes then Coulomb distorted nuclear phase 
shift. In this case the free solutions / and g must be replaced by the appropriate Coulomb wave functions Fi{r], p) and 

(77, p), respectively, and the same must be repeated for the derivatives. The Coulomb corrected scattering length is 
given by a more complicated expression 

A - -{i?$[,(o, i?) - [1 + /:(i?)]$o(o, i?)}/{i?e[,(o, i?) - [1 + /:(i?)]eo(o, R)} (wo) 

where prime denotes the derivative with respect to r and the calculated C{R) takes also account of the Coulomb 
potential. 



B. Bound states 



In the case of a bound state the momentum is imaginary p — > Ik and the wave function will be a function of k and r. 
Our goal is to calculate k which determines the binding energy. In order to enforce the proper threshold behavior of 
the wave function, we set '4'{n, r) = r^'^^ 0(k, r) where the hitherto unknown function ^(k, r), satisfies the Schrodinger 
equation 

d20(K,r)/dr2 = -[2(^+ l)/r] d(/)(K,r)/dr+ + 2/iy(r)] (/)(«, r). (141) 

For V{r) = eq. I|141l) has two linearly independent solutions [nr)^^ ii{Kr) and [nr)^^ ki{Kr). The sought for 
function (/)(k, r) is the regular solution of (|141|l decaying exponentially for large r. It will be convenient replacing the 
single second order equation H141|) by a system of two, first order differential equations, which is accomplished by 
introducing formally a second function x('^: ^) ^-nd equal to the derivative of the wave function x{^i = d0(K, r)/dr. 
The wave equation H141|l . takes the matrix form 



(/.'(K,r)\ / 1 \ /0(K,r) 

X'{^,r)) \ [k-" + 2^,V{r)\ -2(^ + l)/rj 



(142) 



Although the variable r varies between zero and infinity, for practical reasons infinity must be replaced by some cutoff 
radius R which is much bigger than the range of the potential. Thus, in the following it will be assumed that r belongs 
to the interval [0, i?]. We are interested in the regular solution of H141() and in view of the fact that the threshold 
factor controls the rate of the fall off down to zero of the wave function, (j>{K,r) must go to a constant for r ~> 0. 
Since this constant would be absorbed in the normalization, we take its value to be equal to unity. At r = we have 
x(k, 0) — c where c = /i lim,.^o rV{r)/{£ + 1). 

Integrating eq. (|142|) in the limits [0, R] and introducing the Chebyshev grid, we are presented with a system of 
algebraic equations 



1 -^RW- \ f[^]\ _ f[i] 

-^RW-o[p] l+^^RW-o[q]) {[x] 



(143) 
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where [p]i = + 2/Lt V([a;]i); [q]i = 2{i+l)/[x]i. When the system of hnear algebraic equations H143() has been solved, 
the wave function and its derivative can be reconstructed at any point from and (|137|l . respectively making the 

substitution p ^ k. 

To locate a bound state we impose the requirement that the wave function falls off exponentially at large separations. 
In the asymptotic region r > R the potential V{r) becomes negligible and the wave function must be a superposition 
of the free solutions 

(^(k, r) w (Kr)-''' [A{k) kiinr) + B{k) ii(Kr)], (144) 

where the two coefficients A{k) and B{k) may be determined by matching at r = i? the form l|144|l and its derivative 
to the functions H13t)|) and H137|) . respectively. This gives two algebraic equations from which A{k) and B{k) are 
calculated. A bound state wave function must be square integrable which implies that the exploding term proportional 
to in l|144() is inadmissible and therefore the bound state condition ensuring exponential fall-off, is B{k) = 0. 
Expressing B{k) in terms of the wave function and its derivative, we obtain 

(t){K, R) k'einR) - hinR) [x(k, R) + (j){K, R) i/R] = (145) 

where prime denotes derivative with respect to R. Solving eq. (|145|l for k, we are in the position to locate the bound 
states. 

When the Coulomb interaction is present the asymptotic form of the wave function 1)1 44() needs to be modified, and 
the appropriate formula reads 

(t>{K, r) w (Kr)-^-i [A{k) h{Kr) + B{k) /(/tr)]. (146) 

In the case of a bound state we seek the exponentially decaying solution which implies that the term proportional to 
/ must be absent and this provides the condition for the existence of a bound state 

0(k, R)\£ + 1~kR h'{KR)/h{KR) \ + R x{k, R) = 0. (147) 

The presence of negative energy Coulomb wave functions in (|147ll is only a minor complication owing to the fact that 
all that is needed is the logarithmic derivative of the function h given by H12U|I . The latter combination can be very 
efficiently computed using a continued fraction representation. 



VI. SOLUTION OF THE LIPPMANN-SCHWINGER EQUATION 

The L-S equation is an attractive alternative to the Schrodinger equation. The method of solution of the L- 
S equation in the configuration space presented in this section will be applied exclusively to the Volterra integral 
equation of the second kind. This is true for the continuous and discrete spectrum alike. 



A. Continuous spectrum 

To calculate the phase shift from H109() . all that is needed is the wave function u{p,r) on the Chebyshev mesh 
in the interval [0, i?]. The Fourier-Chebyshev interpolation is never needed because the integrals entering 1)1091) are 
computed by applying the Chebyshev quadrature H7UI) which requires the wave function only at the collocation points. 
The wave function is obtained by solving the Volterra equation (|107|l following closely the algorithm specified in H8U|I. 
Introducing the vector [u] containing the values of u{p, Xi), we are led to a system of algebraic equations 

{l-^R]CoW-}-[u] = [fl 

where the kernel matrix K,, is /C^ = (2/i/p) {fi gj — gi fj) V{xj) and the vectors [/] and [g] represent the free solutions 
f(pxi) and g(jpXi), respectively. 

In some applications, to reduce the storage size, it might be necessary to introduce a composite integration by 
chopping the integration interval [0, R\ into a number of smaller partitions. In order to adapt the above scheme to 
such a situation we shall need a second linearly independent solution of the wave equation that necessarily would be 
irregular at the origin. This solution, hereafter denoted as w{p,r), is defined by the appropriate boundary condition 
at infinity. As a convenient choice we consider the solution of the equation 

w(p, r) = g{pr) H / [f{pr') g{pr) - g{pr') f{pr)] V{r') w{p, r') dr'. (148) 
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The Wronskian W[w, u] does not vanish and, similarly to the free propagation case, is equal to p. Therefore, any 
solution of the wave equation can be written as a superposition of the two linearly independent solutions: u and w. 
We wish now to break the interval [0, R] into M partitions [R\-i, R\\ with A = 1,2, - ■■ M where Rq = and Rm = R. 
For greater clarity, we have reserved hereafter a Greek index A to number the partitions and letting Roman indices to 
label the Chebyshev mesh points. It is important to distinguish between these two kinds of labels, one enumerating the 
sub-mesh points (dimension N) within a given partition, and the other enumerating different partitions (dimension 
M). Obviously, the size of the whole mesh is M x N and, accordingly, the collocation points ought to carry labels of 
both types just mentioned. Since the Greek index A numbers the partitions and Roman indices label the Chebyshev 
mesh points, the collocation points are denoted as with A = 1, 2, . . . , M — 1 and i = 1,2, . . . , N , and we have 



X 



X 



i(i?A + Rx-i) + URx - Rx-i) cos[^(i - \)IN]. (149) 



In every partition A = 1, 2, . . . , Af we have two linearly independent local solutions ux and wx satisfying, respectively, 
the appropriate Volterra equations 

ux{p,r) = f{pr) / [f{pr')g{pr) - g{pr') f{pr)]V{r')ux{p,r')dr'; (150) 

and 

2^ [^^ 

wx{p,r) ^ g{pr) + — / [f{pr')g{pr)~g{pr')f{pr)]V{r')wx{p,r')dr'; (151) 
P Jr 

where it is understood that r belongs to the interval [Rx-i, Rx\- The general solution in the latter interval may be 
written as Axux{p,r) + Bxwx{p,r) where Ax and Bx are two constants to be determined. Obviously, the global 
solution constructed from the local solutions is expected to be a continuous and smooth function which implies that 
the different local solutions, together with their derivatives, should be matched at the partition boundaries r = Rx- 
This gives 2{M — 1) equations for 2M coefficients Ax and Bx 

Axux{p,Rx) + Bx'Wx{p,Rx) = Ax+iux+i{p, Rx) + Bx+iwx+i{p, Rx) 
Ax u\{p, Rx) + Bx w'^ip, Rx) = Ax+i <+i(p, Rx) + Bx+i w'^^^{p, Rx), 

with A = 1, 2, ... , M, resulting in two recurrences 

Ax+i=axAx + bxBx; (152) 
Bx+i ^axAx+ Px Bx (153) 

ax = [ux{p,Rx)w'^_f_i{p,Rx) - u'^{p,Rx)wx+i{p,Rx)]/d; 

bx = [wx{p,Rx)w'^_f_i{p,Rx) ~ w')^{p,Rx)wx+i{p,Rx)]/d; 

"A = [u'^{p,Rx)ux+i{p,Rx) - ux{p,Rx)u'x+i{p,Rx)]/d; 

P\ = [w'x{p,Rx)ux+i{p,Rx) - wx{p,Rx)u'^+i{p,Rx)]/d; 

d = ux+i{p,Rx)w'x+i{p,Rx)-u'x+i{p,Rx)wx+i{p,Rx)- 

The above expression involve the derivatives of the two local solutions ma and wx which are obtained by differentiating 
(|15U|I and p51|l . Although in both cases the integrand vanishes for r' = r, but there would be a contribution from 
the derivative of the kernel. One has 

u\{p,r) = f'ipr) - ^ r [f{pr')g'{pr)-g{pr')f{pr)] V{r') ux{p,r') dr' ; 



with 



P 



w'x{p,r) ^ g'ipr) + — I [f{pr')g'{pr)-g{pr')f'{pr)] V{r') wx{p,r') dr' ; 



and 



P 

where prime denotes derivative with respect to r. 

Since the global solution must be regular, the starting values are Ai = 1 and Bi — and this allows to solve 
the recurrence for all the remaining coefficients. The phase shift is obtained by matching the global solution and its 
derivative at r = i? to the asymptotic form A[f{pr) + tan 5 g(j>r)\ . As a result, the phase shift would be obtained from 
l|138|l putting il>{p, R) = Am um{p, R) + Bm wm{p, R)- 
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B. Bound states 



To locate the bound states we need to calculate the Fredholm determinant A(k) from (|125|l . This, in turn, requires 
the solution of of the Volterra equation of the second kind (|126|l . This equation is solved for the wave function u{k, r) 
on the Chebyshev mesh using the algorithm (|80|l . Introducing the vector [u] containing the values of u(K,Xi), we are 
led to a system of algebraic equations 

{l-^RlCoW-}-[u] = [f], 

where the kernel IC, is /C^ — (2/i/K) {fi hj — hi fj) V{xj) and the arrays [/] and [h] represent, respectively, the grid 
values of the free wave functions f{KXi) and h^KXi). With the solution [u] in hand, the Fredholm determinant, is 
given as 

A(k) = H hi V{xi) Ui Wi. 

K ^ — ' 

i=l 

Finally, the bound states would be determined by finding the zeros of A(k). 

If need arises, composite integration algorithm may also be applied, and the procedure is very similar to that used 
in the continuous spectrum case. Thus, the [0,R] range is divided into M partitions i?2, ■ ■ ■ , Rm and in each 
partition A we have two linearly independent wave functions u\ and wx . These wave functions are obtained by solving 
the appropriate Volterra equations 



u\{K,r) = f{Kr) - 



2fi 



f(nr') h{Kr) — h{Kr') f{Kr) V{r') u\{k, r') dr'; 



(154) 



and 



f{Kr')h{nr)-h{Kr')f{Kr) y(r') wa(k, r') dr'; (155) 
R. The wave function would be a superposition of the 



where r is in the interval [Rx-i-, R\\ with Rq ^ and Rm 
two linearly independent local solutions 

u{k, r) = Ax ux{n, r) + Bx wx{n, r). 

The coefficients Ax and Bx can be determined by smoothly matching the different local pieces of u and its derivative 
with respect to r at all partition boundaries. The matching conditions, are 

Axux{n,R\) + Bxw\{k,Rx) = Ax+iux+iiK,Rx) + Bx+iwx+i{k, Rx) 
Ax u'xin, Rx) + Bx w'x{k, Rx) = Ax+i ua+i(«> -^a) + ^a+i Wx+i{k, R\), 
and yield the following recurrence relations for the coefficients 

Ax+i=axAx + bxBx; (156) 
Bx+i ^axAx + I3x Bx (157) 

where 

aA = [ux{n, Rx) w'x+i[k, Rx) - u'x{k, Rx) wx+i{k, Rx)]/ d; 
bx = [wx{k, Rx) w'x+i{k, R\) - w'x{k, Rx) wx+i{k, Rx)]/ d; 
ax = [u'xiK, R\) ua+i(k, -Ra) - ua(k, -Ra) wa+i(«7 Rx)]/d; 

= [w'xiK, Rx) ux+i{k, Rx) - wx{k, R\) u'x+iiK, Rx)]/ d; 
d = ux+iin, Rx) w'x+iiK, Rx) - u'x+i{k,R\) wx+ii^^Rx)- 

The derivatives occurring in the above expressions are obtained by differentiation of the Volterra equations (| 15411 and 
(|155|l . and we get 



u'x{i^,r) = /'(Kr) 



2/i 

K 

w'x{k, r) — h' {nr) H — — 

K 



f{Kr') h'inr) - h{Kr') /'{nr) V{r') ux{n, r') dr'; 



/{nr') h' {nr) — h(Kr') f' (kt) V(r') wxin, r') dr'. 



(158) 
(159) 
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The bound state function must be regular which impHes that Ai = 1 and Bi = are the proper starting values in 
(|15(j|l and (|157|l . The wave function in the asymptotic region r > R, has the form 

u{K,r) A{k) f{Kr) + I3{K)h{Kr) (160) 

where A{k) and B{k) are obtained by matching at r = i? the asymptotic solution H16U|) to the wave function from the 
last partition and subsequently repeating this procedure for the derivatives. For arbitrary value of n the wave function 
(|160|l would not be square integrable owing to the presence of exponentially increasing function /(k, r). However, for 
some particular k values for which A{k) = the exploding term is absent and the wave function shows exponential 
fall-off for large r. For such k we have a bound state. The ultimate bound state condition is obtained by calculating 
A{k) from the matching condition and equating it to zero. This gives 

[Am ua/(k, R) + Bm wm{k, R)] h'inR) - [Am u\,[{k, R) + Bm Wm{h, R)] h{KR) ^ 0. (161) 

The bound state algorithm is implemented as follows. Firstly, for each of the M partitions the two sets of TV algebraic 
equations must be solved 

{1 - i(i?A - Rx-i) IC^oW-}- Ka] = [fx], 
{1 + i(i?A - Rx-i) ICx o W+} ■ [wx] = [hx]. 
where for fixed A the kernel ICx is a N x N matrix 

[ICxh - i2fi/K) [f{Kxf)h{KX^) - h(Kx':)j(nx])] Vi^x]) 

and the arrays [fx] and [hx] represent, respectively, the grid values of the free wave functions f{Kxf) and h{Kx^) 
with i,j = 1,2, . . . , N. The solutions at the collocation points supply enough information to reconstruct the wave 
functions locally by Chebyshev interpolation, and in particular at the matching points Rx-i and Rx- 

Knowing the local functions ux{K,r) and wx{K,r) at the collocation points that belong to the partition A, we are 
also in the position to compute the r-derivatives of these functions with the aid of (|158|) and H159II by performing the 
Chebyshev integrations. The derivative arrays, are 

[u'x]^[fx]-}C'xoW+-[ux] 
[w'x]^[hx]+IC'xoW+-[wx] 

with the derivative kernel given as 

[/C^]., = i2t,/K) [f{Kxl)h'{KX^) - h{^X^)f'{^xj)] V{xj) 

where prime denotes the derivative with respect to r. 



VII. MOMENTUM SPACE SOLUTIONS 



A. Continuous spectrum 

In order to solve 1)129(1 . we take a simple rational transformation 

q = pcr{l + t)/{l-t) (162) 

mapping the infinite integration range onto the [—1, 1] interval. The presence of an adjustable slope parameter a 
in 1(162(1 improves the flexibility of this transformation. Now, invoking 1(25(1 . both, the potential and the off-shell K- 
matrix, will be expressed in terms of the Chebyshev cardinal functions. For the potential we just use the interpolating 
formula 

N 

Ut{k', k") = G,{t') G,{t") (163) 
with Uij = Ui{ki, kj) and kj ~ pcr(l + tj)/{l — tj). Similar expansion is used for the off-shell K-matrix 

N 

{k'[Ke{p')[k") = J2 G,{t')K,,{p')G,{n (164) 
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where in both cases above the rational mapping (|162|l provides the relation connecting k' with t' and k" with t", 
respectively. The hitherto unknown matrix Kij(j>^) will be determined by solving H129|) . Using the new variables and 
the expansions H163|l and p64(l . the L-S equation takes the matrix form 

(l + U-T)- K{p^) = U (165) 

where the matrix F is given as the integral 

with T = (1 — (t)/(1 + tr). The integral occurring in Ij 166(1 will be evaluated using the automatic quadrature ((47(1 . The 
presence of the cardinal functions in the integrand brings major simplifications because in the summation over the 
Gauss-Chebyshev abscissas only a single term survives and the resulting matrix T is diagonal 

r,y = 5y 4paV[^(l + w,(r)/[l - Ut] [(1 + t,)/(l - U)f . 

The explicit solution of p65(l 

Kip"^) = (1 + u -ry^ -u (167) 

inserted in 1164(1 . yields the desired off-shell K-matrix. 

It is evident from ((162(1 that the mapping introduces a singularity at t = 1 reflecting the fact that the original 
integral was in the infinite range. The singular point t — 1 corresponds to q going to infinity in ((1291) but, in all cases 
of physical interest, the potential goes to zero much faster than quadratically when this limit is taken, making the 
t = 1 singularity quite harmless. 

When only the phase shift is needed it is sufficient to calculate the half-off-shell K-matrix {kj \Ki{p'^)\p) on the 
Chebyshev mesh kj, j — 1,2, ... N. The on-shell K-matrix is subsequently obtained, from the expansion 

N 

(p|if,(p2)|p)=^G,(r) {k,\Ke{p')\p), 

where 

N 

{k, mp')\p) = ^ (1 1/ • r)-/ u{h,p). 

The linear system ((165(1 needs to be solved only once to get the on-shell K-matrix whereas to determine the fully 
off-shell K-matrix this operation must be repeated N times. 



B. Bound states 



In momentum space the Schrodringer and the L-S equation both result in the same scheme for locating the bound 
states. We begin with the Schrodinger equation and we wish to apply the Gauss-Chebyshev quadrature which requires 
that the [0,oo] integration range to be mapped onto [—1, 1] interval. This may be accomplished by a change of the 
integration variable 

k' = {a/a){l + t)/{l-i) (168) 

where a is the range of the potential and a is an adjustable slope parameter. The integration over t is then carried 
out by the automatic quadrature 1(43(1 and we end up with a real and symmetric eigenvalue problem 

[1 - \{k) M{k^)] ■ u = Q (169) 

which can be solved by standard methods. The explicit form of the matrix M., is 



M^j{k^) = -i^a/a^) [V^/(l - U)] {hJ^K^ + k1) Ut{h,k,) {k,l^^^]) - t,)\ 
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where the new wave function is uj = [y/Wj/{l — tj)] u{kj). Solving the algebraic eigenvalue problem repeatedly for a 
number of k values, we get in each case the value of X{k). When the latter function has been tabulated, it may be 
locally approximated by polynomials and it remains to apply inverse interpolation to solve the equation A(ko) — 1- 

The above procedure applies also for Coulomb potential but this case requires some care. The point charge Coulomb 
potential in momentum space Ui{k, k') exhibits in all partial waves a logarithmic singularity aX k' — k (cf. Appendix). 
As a result, the diagonal elements of the M. matrix become infinite and the momentum space framework is in trouble. 
In the literature, this difficulty is alleviated by a sophisticated subtraction scheme serving to eliminate the singularity 
j34j |35l| . Within the Chebyshev approach, however, the Coulomb problem in momentum space can be solved directly 
without subtractions. In order to achieve this goal the integrals involving the logarithmic singularity need to be 
approximated by the dedicated quadrature H57I) . 

It might be instructive to dwell on the momentum space hydrogen-like problem. In this case, we are presented with 
the integral equation (for the notation cf. Appendix) 



2 

u{0 = - / {Pi{z) log(e' +0 - Wi.,{z) - P,{z) logic' -^1} 

71" Jo 



w(f')de' 



(170) 



In (jl7()|l we switched to dimensionless variables, and C, and x are defined, respectively, as ka, k'a and na where a 
is the appropriate Bohr radius. The integrals in 1)170(1 containing the non-singular terms can be approximated by the 
Gauss-Chebyshev quadrature, whereas the last integral involving the singular log |C' — ^\ term requires the dedicated 
quadrature H57() . Changing the integration variable 

C = <j(l+t)/{l-t) (171) 

the resulting matrix M , takes the form 

. ^ 4ct wj {Pejzjj) log |1 - tj tj\ - Wi^i{z,j)} - PeXz.j) »j(t^) ^^^^^ 

where = (£,f + / and it is apparent that the diagonal terms are all finite and well defined. Since the 
matrix M in (|172|l is not symmetric anymore the eigenvalues have to be determined by solving the secular equation 
det |1 — A1(a;^)| = 0. The same procedure as above applies also to other problems in which the point charge Coulomb 
potential is supplemented by an additional short range potential. 

We turn now to the approach based on the L-S equation. When the momentum p is on the imaginary axis, we set 
as before p = in and the L-S equation for the T-matrix, reads 

{k'\TAn')\k)=Ue{k',k)-- rUe{k',q)^^{q\Te{K')\k) (173) 

and the integral equation is non-singular. We are going to use again the same mapping (|168|l . Expanding the potential 
and the T-matrix in the Chebyshev cardinal functions, we set 

N 

{k'\Te{K')\k") = G,{t')T,,{^)G,{t") (174) 

and the L-S equation H173|) takes the matrix form 

T(k2) = U -U- r(K2) •T(k2) (175) 

where the matrix r(K^), is 

2^ 4cr Gi{t)Gj{t) dt 



^ " i_i 1 + {na/af{l - t)2/(l + t)^ (1 - ' ^^^^^ 
Applying Gauss-Chebyshev quadrature, we obtain 

r.y(K') = S,j (Aa/na) w,/{l - t,)V[l + (W^)'(l - i,)V(l (177) 

1 1 

and it is evident that the matrix T is positive definite. Therefore, we may dissolve F as Fa • r2 in ((175(1 and defining 
1 1 

a new t-matrix T = F2 • T • F2 , the L-S equation, becomes 

{1-M(k^)} ■T{k^)^~M{k^), (178) 
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where M = —Fa • U • T^. Since the calculation of the t-matrix involves inverting the matrix 1 — A^(k^), the poles of 
T will occur at such values of k at which the determinant of the latter matrix vanishes. Thus, the ultimate equation 
for K, takes the form 

Det|l - = 0. (179) 

Clearly, (|179|l is nothing else but the secular equation associated with the eigenvalue problem Hlt)9|) which demonstrates 
the equivalence of L-S and Schrodinger approach. 

The above scheme can be also used to calculate the scattering length A. From (|174|l . we obtain 

N 

A = -{0\T,{0)\0) = G.(-1)T,,(0)G,(-1) (180) 

with 

T{0) = [1 + U -TiO)]-^ -U. (181) 

The matrix r(0) is diagonal and takes a simple form 

r,,(0) = % (4s/7r) [aw, /{I - t,f]. (182) 

Summarizing, the scattering length is given as a scalar product of two vectors 

A==-[G]-[X] (183) 

where the cardinal functions in the array [G] are to be taken at i = — 1 and the vector [X] is obtained as the solution 
of a linear system of equations 

{1 + U -riO)} ■ [X] = U ■ [G]. (184) 



VIII. NUMERICAL PERFORMANCE TEST 



Having set up the theoretical background, we are ready now to investigate the performance of the semi-spectral 
Chebyshev method by solving numerically some concrete differential and integral equations encountered in quantum 
mechanics and subsequently comparing the results with the exact solutions. As mentioned before, we confine our 
attention to time independent non-relativistic single-channel two-body problems, disregarding all internal degrees of 
freedom. We considered three popular potential shapes: (i) exponential, (ii) Hulthen and (iii) Morse for which exact 
solutions are available and have been collected in the Appendix. We solved numerically the s-wave scattering problem 
and located the bound states for these potentials by solving three different equations: (i) the differential equation 
(Schrodinger), (ii) the Volterra integral equation (L-S in configuration space), and (iii) the Fredholm integral equation 
(L-S in momentum space). Finally, we solved numerically the bound state problem for a point-charge Coulomb 
potential, both in the configuration and in the momentum space. 

We begin with the continuous spectrum where we shall calculate the s-wave phase-shift for three different potentials 
for which exact expressions are available (cf. Appendix). As a result of scaling, the phase-shift becomes a function 
of two dimensionless variables 5{s,S.) where ^ is the center-of-mass momentum in inverse range units (^ — pa) and s 
is the potential strength parameter. Since in each case we do know the exact phase-shift, the quantity of interest for 
us here will be the relative error. Our intention was to remain close to the realm of nuclear physics and the strength 
parameter s has been chosen in such a way that the potentials roughly reproduce the proton- proton ^5*0 phase-shift. 
This interaction is quite strong, being almost capable of supporting a bound state. The dependence upon s will be 
examined later on. For each of the considered here potentials, keeping s fixed, we have computed the phase-shifts at 
Up equidistant values of the momentum i = 1, 2, • • • rip and this has been done repeatedly for gradually increasing 
approximation order N. For an assigned approximation order N, we may define an average relative error E(N), given 
as 

where S and 5n denote, respectively, the exact and the N-th approximant to the computed phase-shift. In the actual 
computations we took n„ = 100 with s = 0.8 for both, exponential and Hulthen potential and s = 0.2 for Morse 
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FIG. 1: Phase-shift (L=0) for different potentials FIG. 2: Average relative error H185|l for the phase- 

vs. momentum in inverse range units. shift versus A'' (exponential potential). 



potential. The corresponding phase shifts vs. the center-of-mass momentum ^ are presented in fig. ^ Among the 
considered potentials, only the Morse potential, which is repulsive at small separations, is capable of producing a 
cross-over in the phase-shift. The error defined in l|185|) depends upon the order of the semi-spectral approximation 
A'^ used in the computation of (5 at and the quality of the approximation will be determined by the rate at which the 
error goes to zero for large TV. The function E{N) reflecting the convergence of the semi-spectral method is displayed 
in figs. 13- El for the three different potentials. 

In the procedure described above the value of the potential strength was frozen while we were calculating the 
phase-shifts at the different momenta. Now, we wish to examine the opposite situation when the momentum is fixed 
and the strength is allowed to vary. However, instead of the phase-shift we shall consider the scattering length which 
is known to posses a complicated structure as a function of s. This function exhibits interlacing zeros and poles and 
changes by many order of magnitude for s values close to a pole, assuming both, positive and negative values. The 
plots of the scattering length vs. s for different potentials are presented in fig. |31 Similarly as before, we consider Up 
values of the strength s^; i = 1, 2, ... Tip and introduce the error function 

""^^^^^g — — ' ^'''^ 

where A{si) and AN{si) denote, respectively, the exact scattering length (cf. Appendix) and its N-th approximant 
computed by the semi-spectral method. In (|186|l we just take the arithmetic mean of the relative errors computed 
for different s. In the actual computations we took Up — 100 adopting for Si an equidistant sequence of values in the 
intervals (0, 2.5), (0, 1.1), (0, 0.3) for the exponential, Hulten and Morse potential, respectively. As seen from fig. 
|5l in each case the considered interval contained the s value at which the scattering length exhibits a pole. These 
critical values of the strength are s « 1.44577 for the exponential, s = 1 for the Hulthen and s = i for the Morse 
potential, respectively. The error function obtained from H186(l for different potentials is given in figs. I6I8I In 
the discrete spectrum problem the dimensionless quantity of interest is the imaginary part of the momentum (in 
inverse range units) at which the T-matrix has a pole. This, necessarily non-negative quantity, will be denoted as 
x{s) where we have emphasised its dependence upon the potential strength s. Thus, similarly as for the scattering 
length A(s), a function of s needs to be determined but it is much harder to obtain the exact solution. The function 
x{s) is known explicitly for the Hulthen potential but (cf. Appendix) for the exponential and the Morse potential 
x{s) is provided in an entangled form, as a solution of a transcendental equation of the form F{x, s) = where 
the function F might be quite complicated. In general, such equation can be solved only numerically and this is 
a difficult task especially when machine accuracy is desired. Nevertheless, for the exponential potential for some 
particular x values, s{x) can be obtained quite easily and one such obvious instance is a; = 4 in which case we get 
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FIG. 3: Average relative error (11851 for the phase- FIG. 4: Average relative error 118511 for the phase- 

shift versus A'^ (Hulthen potential). shift versus N (Morse potential). 




FIG. 5: Scattering length in inverse range units vs. 
s for different potentials. 



FIG. 6; Average relative error 11861 for the scat- 
tering length versus A'' (exponential potential). 
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FIG. 7: Average relative error 118611 for the scat- 
tering length versus A'^ (Hulthen potential). 




FIG. 8; Average relative error 11861 for the scat- 
tering length versus A'^ (Morse potential). 
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FIG. 11: Relative error for the calculated binding 
energy - Morse potential. 



FIG. 12: Relative error for the calculated binding 
energies - Coulomb potential. 



s — (7r/2)^ as the exact solution. Unfortunately, for the Morse potential there is only the hard way left. For testing 
purposes of the semi-spectral method as the benchmark solutions we used the exact values s — (7r/2)^, x — j for 
exponential potential and, respectively, s = 3/2, x — 1 for the Hulthen potential. For the Morse potential the 
benchmark solution was obtained numerically. The shape was adopted from @ where Morse potential was used as a 
model of proton-neutron interaction in the deuteron (with a — 0.3408 fm and d = 0.8668 fm). We slightly retouched 
the depth so that the potential better reproduces the deuteron binding energy B — —2.22 MeV and at the pole, we 
have s = 0.33509414149514, x = 0.0078864302204068. The above pairs of values of {s,x) for different potentials 
were regarded as exact for testing purposes of the semi-spectral method. Using the s value as input, the values of 
X was subsequently determined by solving the Schrodinger equation, and the L-S equation in the coordinate and in 
momentum space, respectively. The relative error for x is presented in Figs. 1911 II as a function of the approximation 
order TV. Unexpectedly, a comparison of the results obtained by different methods in Figs. I2I11I shows that the 
unpopular Volterra equation method appears to be the winner. 

Finally, in Fig. 1121 we display the relative errors for the calculated binding energies in a hydrogen-like system where 
a point-like charge Coulomb potential was operative. The values of the quantum numbers (n, £) are given on the 
plot. We wish to note that in the momentum space calculation there was no need to invoke Lande's subtraction 
technique and the difficulty connected with the Coulomb singularity was avoided by using the dedicated weights 
r2j (z), as explained in detail in Sec. VII. Therefore, the results presented in Fig. 1121 indirectly involve also a test on 
the accuracy of the method applied for calculating the singular integrals. 

We have been displaying the logarithm of the error versus N as this is convenient for a quick assessment of the 
asymptotic convergence rate and a linear plot would indicate a geometric convergence (as a reminder we note that 
an algebraic convergence would lead to a linear dependence between log(£') and log(A^)). It has to be kept in mind, 
however, that with N < 100, as adopted in our computations, the approximation order N cannot be really regarded 
as asymptotic {N — > oo) and a more complicated dependence upon N should not be surprising. Nevertheless, as 
might have been expected, the semi-spectral errors presented in this Section decrease very rapidly to zero, in most 
cases showing indeed a geometric convergence. 



IX. SUMMARY AND CONCLUSIONS 



The purpose of this work was a concise presentation of the mathematical techniques associated with the Chebyshev 
pseudo-spectral method. The emphasis was on the practical aspects in the efficient implementation of this method on 
high-performance computers, in a manner accessible to newcomers to the field. The paper is self-contained, i.e. apart 
from a standard linear algebra package, no other resources are needed and the supplied algorithms can be immediately 
turned into a suite of computer codes capable of solving a wide variety of practical problems. Formal considerations 
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have been illustrated by providing concrete solutions of some typical quantum mechanical problems. 

In our experience the pseudo-spectral Chebyshev method seems to be without downsides: the implementation and 
programing could not have been simpler, the precision and stability are excellent with the convergence being usually 
exponential. Although the resulting matrices are not sparse but it does not matter since their sizes are not large. 
Finally, the same scheme can be used for solving differential, integral or integro-differential equations. 



APPENDIX A: CONFIGURATION SPACE 

Denoting the range as a and the depth as Vq, we introduce two dimensionless quantities: the potential strength 
parameter s = 2/i Vq and the momentum in inverse range units ^ = pa. The s-wave phase shift S{s, ^) is a function 
of s and ^ and the scattering length in the units of range A{s)/a is defined as the limit of S{s,£^)/^ when ^ 0. 

(a) Exponential potential. 

2nV{r) = -(s/a2)e-''/« 
and only the attractive case will be considered s > 0. The phase shift, is 

S{s, In [ J2j^(2v^) r(l + 21^)1 - e In s 

where Jy{x) denotes a complex order v Bessel function of the first kind. The scattering length, is 

A(s)/a= -2(7 + logV^) + ^ro(2V^)/Jo(2yi), 

where 7 is Euler constant and Yq{x) is the zero-order Bessel function of the second kind. For an assigned value of 
s, the bound states poles are at ^ = ix and are obtained by solving the transcendental equation 



(b) Hulthen potential. 



■V{r) ^-{s/a^) (e^/«-iy 



and only the attractive case will be considered s > 0. The phase shift, is 



-3 inr(i + 2iO + inr(v/I^-iO + inr(-v^^- 

and the scattering length, is 

Ais)/a = 27 - V^l + \/^) - ^{1 - V^), 
where ip denotes the digamma function. The bound state poles occur at ^ = ix where 

a; = (s - n2)/(2n), n= 1,2,3... 

(c) Morse potential and Morse barrier. 

In both cases the potential is given by the same expression 

Vir) = -(s/a^) e(^ - '^)/« [2 - e(^ " '^)/« 

where d is a parameter of the dimension of length. For Morse potential we have s > 0, whereas s < for Morse 
barrier. The phase shift, is 



argM(i + ie- V^, l + 2iO z), 



s > 
s < 



arg M(i + iC - iV^, 1 + 2iC, iz)-^z, 

where z = 2e'^/^ \/W\ ^-nd M(a, 6, z) denotes the Kummer function 4]. The scattering length, takes the form 
A{s)/a = -27 - V(i - Vi) - log z - r(i - Vi) C/(i - Vi, 1, z)/Af (i - ^/i, 1, z) 
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for s > and 

A{s)/a = -27- SRjV'lj -logiz + r(i - i^/^) U{\ 1, iz) / M {\ - , 1, \z)] 

for s < where C/(a, 6, z) denotes the Tricomi hypergeometric function 4J. For assigned values of s and d/a, the 
bound state poles are at ^ = ia; where x is a solutions of the equation 

M(i +x- Vs, 1 + 2a:;, z) = 0, s>0 
5RM(i + a; - i^/^, 1 + 2a;, iz) = 0, s < 0. 

(d) Coulomb potential. 

For a Coulomb potential of the form V{r) = —a/r where a is the fine structure constant, we have 

2^JiV{r)^-{2/a^){a/r) 

where a = l/iia is the Bohr radius and only the attractive case will be considered. The phase shift, is 

5,(0 = argr(£ + l-i/C). 

Bound states occur at ^ = i a; where 

x = l/{n + £+l), n 0, 1,2, . . . . 
When a Z factor multiplies the Coulomb potential, as result of scaling r r/Z we get x ^ Z x. 

APPENDIX B: MOMENTUM SPACE 

The potentials in momentum space are obtained from H128|) . Introducing the abbreviations x = a{k + k'), y = 
a{k — k'), we have 

(a) Exponential potential: 

Uo{k,k')^-2sa/[{l + x^){l + y% 

(b) Hulthen potential: 

Uoik, k') = -2sa ^^i^ + ^^l-mi + iy) ^ £ , 



(c) Morse potential and Morse barrier: 



Uo{k, k') = —sa 



(l + a;2)(l + y2) ' ^ 4 [1 + ( 1 ^)2] + ( 1 j^)2] ' 
(d) Coulomb potential: 

Ue{k,k') = ^4a Q,(z)/(a;2 - y2) 

where z = (a;^ + y'^)/{x'^ — y^). Following Qi{z) denotes the Legendre function of the second kind which 
exhibits a logarithmic singularity 

Q,{z)^P,{z) ilog^-M/,_i(z) 

with W-i{z) = and 

n 



Wi-i{z) = V ip„_i(z)P^_„(z), 
where Pi{z) denotes the Legendre polynomial. 
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